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1.  INTRODUCTION 


The  term  "blind  deconvolution”  refers  to  methods  that  jointly  estimate  an  object  and  a  blurring  function 
from  a  measurement  that  consists  of  their  convolution  (usually  corrupted  by  noise  as  well).  It  is 
straightforward  to  show  that  the  blind  deconvolution  problem  does  not  have  a  unique  solution  without 
including  additional  information  in  the  estimation  process.  The  additional  information  that  can  or  should 
be  included  depends  on  the  physics  of  the  measurement  and  the  available  a  priori  knowledge  about  the 
object  and  the  blurring  function.  There  exists  a  plethora  of  blind  deconvolution  publications  describing 
algorithms  and/or  theoretical  results,  differing  primarily  in  the  additional  information  assumed  to  be 
known  and  the  means  by  which  the  object  and  the  blurring  function  are  estimated.  A  review  article  and 
its  companion,  published  in  1996,  provide  a  nice  taxonomy  of  the  types  of  additional  knowledge 
commonly  used  and  the  algorithms  employed  to  carry  out  the  estimation  process  [1,2].  Although  many 
years  have  passed  since  the  articles  were  published,  they  are  still  excellent  reviews  for  those  unfamiliar 
with  blind  deconvolution. 

The  use  of  the  term  "blind  deconvolution”  occurred  first  in  Ref.  3  to  describe  methods  for 
restoring  images  and  audio  recordings.  These  methods  require  a  separate  set  of  measurements  of  so- 
called  prototype  objects  that  share  statistical  similarities  with  the  object  of  interest  in  the  measurements  in 
order  to  recover  both  the  blurring  function  and  the  object.  Although  these  methods  generalized  the 
standard  deconvolution  problem  (where  the  blurring  function  is  assumed  to  be  known),  the  requirement  to 
have  a  set  of  prototype  images  limited  the  applicability  of  the  methods.  A  number  of  years  later,  it  was 
shown  that  the  blind  deconvolution  problem  for  signals  of  dimension  greater  than  one  is  uniquely 
solvable  (up  to  a  scale  factor  and  a  spatial  shift)  using  only  the  information  that  the  object  and  the  blurring 
function  have  finite  supports  [4].  This  result  eliminated  the  need  for  a  separate  set  of  prototype  images. 

One  important  application  of  finite-support-based  blind  deconvolution  is  high-resolution 
astronomical  imaging  using  ground-based  telescopes,  for  at  least  two  reasons.  The  first  reason  is  that 
astronomical  objects  have  finite  size  and  are  on  a  black  background,  and  thus  have  finite  support.  The 
second  reason  is  that  the  structure  of  the  blurring  functions  in  short-exposure  atmospherically-blurred 
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images  is  such  that  the  deconvolution  operation  can  bring  about  significant  increases  in  resolution.  The 
blurring  function,  or  point  spread  function  (PSF),  is  dominated  by  atmospheric  turbulence  effects. 
Atmospheric  turbulence  blurring  is  the  result  of  spatial  and  temporal  fluctuations  in  the  index  of 
refraction  of  the  atmosphere.  In  long-exposure  images  (integration  times  of  seconds  or  longer),  multiple 
atmospheric  turbulence  realizations  combine  to  produce  a  long-exposure  PSF  that  is  essentially 
bandlimited  to  angular  frequencies  less  than  r^/X.  where  is  the  observation  wavelength  and  r„  is  the 
Fried  parameter  [5]  that  has  a  nominal  value  of  10-20  cm  for  visible  wavelengths  at  astronomical 
telescope  sites.  Because  of  the  structure  of  this  PSF.  deconvolution  methods  produce  only  nominal 
improvements  in  resolution.  Furthermore,  since  the  diffraction-limited  resolution  of  telescopes  is  D/X. 
where  D  is  the  diameter  of  the  primary  mirror  of  the  telescope,  and  since  modem  astronomical  telescopes 
have  diameters  of  8  meters  or  more,  long-exposure  PSFs  dramatically  reduce  the  theoretical  resolution 
limits  of  telescopes. 

The  situation  is  significantly  different  for  short-exposure  PSFs  (where  the  atmospheric  turbulence 
remains  fixed  during  the  exposure  time)  because  they  have  a  structure  that  permits  deconvolution-based 
resolution  increases  of  up  to  a  factor  of  D/ro  more  than  is  present  in  long-exposure  images.  The  first 
published  deconvolution  method  to  exploit  this  structure  of  short-exposure  PSFs  (Labeyrie’s  technique) 
generated  deconvolved  energy  spectra  with  diffraction-limited  resolution  [6].  Labeyrie's  technique  was 
then  extended  to  other  techniques  (such  as  cross-spectrum  [7]  and  bispectrum  [8]  estimation)  that  permit 
the  recovery  of  Fourier  phase  spectra  as  well  as  energy  spectra  so  that  diffraction-limited  images  can  be 
obtained  [9].  All  of  these  techniques  (called  speckle  imaging  techniques)  require  as  input  a  set  of  short- 
exposure  images  of  an  object  where  the  object  is  common  to  all  the  short  exposures  but  the  atmospheric 
blurring  is  different  for  each  image.  They  then  estimate  average  values  of  the  above-mentioned  quantities 
for  which  diffraction-limited  information  is  retained  in  the  averaging  process.  The  object  Fourier  phase 
can  be  extracted  from  the  average  image  cross  spectrum  or  bispectrum  without  deconvolution,  but  the 
image  energy  spectrum  is  equal  to  the  object  energy  spectrum  multiplied  by  the  atmosphere  energy 
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spectrum.  Thus  it  is  necessary  to  collect  a  separate  set  of  measurements  of  an  unresolved  star  to  generate 
an  atmosphere  energy  spectrum  estimate  to  be  used  for  deconvolution.  For  this  reason,  speckle  imaging 
techniques  have  much  in  common  with  the  techniques  first  referred  to  as  blind  deconvolution  techniques 
[3],  where  the  separate  star  measurements  are  the  equivalent  of  the  prototype  image  measurements. 

In  the  field  of  astronomical  imaging,  the  term  "blind  deconvolution”  refers  to  deconvolution 
techniques  that  do  not  require  separate  measurements  to  determine  the  atmospheric  PSF.  The  first  blind 
deconvolution  algorithm  applied  to  atmospherically-blurred  images,  to  the  authors’  knowledge,  was  the 
seminal  work  of  Ayers  and  Dainty  [10],  where  they  used  an  iterative  implementation  reminiscent  of  the 
Gerchber  /Saxton  phase  retrieval  algorithm  [1 1]  and  enforced  a  positivity  constraint  as  well  as  a  support 
constraint.  The  Ayers/Dainty  algorithm  was  quickly  followed  by  a  version  that  works  with  complex  as 
well  as  real  images  [12].  Because  the  convergence  and  noise  properties  of  the  Ayers/Dainty  algorithm  are 
not  optimal,  many  other  blind  deconvolution  implementations  for  removing  atmospheric  turbulence 
blurring  from  image  data  alone  have  been  proposed,  including  projections  on  convex  sets  [13],  neural 
networks  [14],  maximum-likelihood  estimation  (MLE)  [15],  non-negative  matrix  factorization  [16], 
simulated  annealing  [17],  and  maximum  a  prior  estimation  [18].  It  was  soon  realized  that  additional 
datatypes  that  were  collected  simultaneously  with  the  image-plane  data  can  produce  higher-quality 
restorations.  Such  datatypes  include  a  separate  image-plane  measurements  that  have  purposely- induced 
phase  distortions  such  as  defocus  [19],  wavefront-sensor  data  [18],  and  imaging  at  multiple  wavelengths 
[20]. 

A  particularly-fruitful  approach  to  blind  deconvolution  is  to  use  a  set  of  short-exposure  images 
(measurement  frames)  of  an  object  where  the  object  is  common  to  all  the  images  but  the  blurring  is 
different  for  each  image.  This  type  of  blind  deconvolution  is  commonly  called  multi-frame  blind 
deconvolution  (MFBD).  The  usual  way  to  generate  this  set  of  images  is  to  collect  a  temporal  series  of 
short-exposure  images  of  the  object.  The  first  use  of  the  MFBD  term  was  in  the  context  of  jointly 
estimating  both  the  unblurred  object  and  the  PSFs  for  each  of  the  measurement  frames  [21],  but  this  was 
not  the  first  approach  to  using  this  type  of  data.  The  first  approach  was  in  the  context  of  speckle  imaging. 
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where  the  energy  spectrum  of  the  system  was  calculated  from  the  measurement  frames  instead  of  using  a 
separate  set  of  measurements  of  an  unresolved  star  [22].  Shortly  thereafter,  an  approach  was  published 
that  used  a  set  of  short-e.xposure  images  sequentially  in  a  modified  Ayers/Dainty  algorithm  [23],  and  then 
extended  to  use  the  short-exposure  images  jointly  [24].  A  key  benefit  of  the  approach  described  in  Ref. 
21  is  that  a  penalized  maximum-likelihood  algorithm  is  used,  which  tends  to  be  more  robust  than 
Gerchber/Saxton  algorithms  or  projections  onto  convex  sets. 

The  vast  majority  of  the  research  in  blind  deconvolution  is  concentrated  in  four  different  areas: 
(1)  algorithm  development,  (2)  types  of  additional  knowledge  to  include  in  algorithms,  (3)  theoretical 
results  on  speed  and  convergence,  and  (4)  applications  to  real  data.  In  addition  to  these  areas,  another 
important  area  of  blind  deconvolution  research  is  the  development  of  algorithm-independent  image- 
quality  performance  bounds  such  as  Cramer-Rao  lower  bounds  (CRBs)  [25].  A  number  of  papers  have 
been  published  including  the  use  of  CRBs  to  compare  two  algorithms  assuming  white  Gaussian  noise  and 
single-input  single-output  systems  [26],  asymptotic  CRB  expressions  for  a  multiple-input  multiple-output 
system  with  white  non-Gaussian  noise  [27,28]  and  analyses  of  phase-diversity-based  blind  deconvolution 
[29,30]. 

A  second  important  area  of  blind  deconvolution  research  for  which  few  published  results  are 
available  is  investigating  algorithm  architectures  that  take  advantage  of  parallel  computational  resources 
and  demonstrations  of  improved  performance.  One  published  approach  is  based  on  segmenting  the 
blurred  images  into  a  number  of  smaller  images,  each  of  which  can  be  processed  in  parallel  [14].  The 
authors  do  not  discuss  the  speed  improvements  achieved  by  this  type  of  parallelization  in  any  detail.  A 
second  approach  to  blind  and  parallel  processing  is  to  set  up  the  conjugate  gradient  search  algorithm  and 
the  filter  bank  structure  to  be  in  parallel  [31],  Again,  no  speed  improvements  are  discussed.  A  third 
approach  is  based  upon  specialized  hardware;  however,  no  discussion  of  how  to  implement  blind 
deconvolution  on  this  architecture  is  given  [32]. 

In  this  paper,  we  describe  our  MFBD  algorithm  that  we  call  the  physically-constrained  iterative 
deconvolution  (PCID)  algorithm.  It  is  a  MLE-based  algorithm  that  utilizes  a  conjugate-gradient 
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minimization  routine  and  regularization  with  the  use  of  either  a  penalty  term  or  a  Fourier-domain  filter. 
Because  our  application  is  astronomical  imaging,  we  have  little  or  no  a  priori  information  about  the 
objects  to  be  imaged  other  than  generic  constraints  such  as  support  and  positivity  constraints;  therefore, 
the  object  is  always  estimated  pixel-by-pixel  in  the  image  domain.  Because  we  do  assume  that  the  system 
PSF  results  from  atmospheric  turbulence  and  telescope  diffraction,  we  have  included  the  ability  to 
estimate  the  PSFs  either  pixel-by-pixel  in  the  image  domain  or  in  terms  of  a  Zernike-based  expansion  of 
the  phase  in  the  pupil  of  the  telescope  [9).  We  also  describe  how  we  parallelized  PCID  to  run  on 
distributed-memory  parallel  computers  in  a  manner  that  scales  well  for  processor  counts  of  one  hundred 
or  more,  depending  upon  the  speed  of  the  communication  network  between  the  nodes,  the  number  of 
pixels  in  a  measurement  frame,  and  the  number  of  measurement  frames  used  in  a  single  image  restoration. 
We  show  that  image  restorations  can  be  generated  in  seconds  to  minutes  for  using  a  Cray  XDI 
supercomputer  [33].  We  demonstrate  that  the  variances  of  PC  I  D-generated  object  estimates  achieve  or 
closely  approach  their  CRBs  for  a  variety  of  algorithm  parameters. 

The  outline  of  the  paper  is  as  follows:  in  Section  2  we  describe  all  but  the  parallelized 
architecture  of  the  PCID  algorithm,  in  Section  3  we  describe  the  parallel  architecture  of  PCID  and 
demonstrate  PCID’s  performance  as  a  function  of  the  number  of  processors  used  by  the  algorithm,  in 
Section  4  we  present  CRB  theory  as  it  applies  to  PCID  and  compare  PCID  performance  to  these  CRBs,  in 
Section  5  we  give  several  examples  of  PCID  restorations,  and  in  Section  6  we  give  conclusions  and 
anticipated  future  work. 

2,  PCID  ALGORITHM  DESCRIPTION 

In  this  section  we  first  present  the  theory  underlying  the  PCID  algorithm  and  then  describe  all  the 
implementation  issues  associated  with  producing  the  desired  estimates.  The  PCID  algorithm  incorporates 
features  that  increase  the  probability  and  speed  of  finding  the  global  minimum,  achieve  superresolution 
by  de-aliasing,  handle  non-idealities  such  as  clipped  images  due  to  jitter  and/or  object  size  relative  to  the 
detector  size,  and  generate  many  of  the  PCID  inputs  automatically. 
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A.  PCID  theory 

The  forward  model  used  in  the  PCID  algorithm  is  the  standard  linear  imaging  model  given  by 

'„W  =  o(x)*//„(x)+«„(x).  (1) 

where  /„,(x)  is  the  m"'  measurement  frame,  o(x)  is  the  true  object,  /i„,(x)  is  the  m"’  PSF,  is  the  m"’ 
noise  realization,  x  is  a  two-dimensional  spatial  location  vector,  bold-face  type  indicates  vector  and 
matrix  quantities,  and  *  denotes  convolution.  The  noise  term  incorporates  both  Gaussian  and  Poisson 
noises  that  are  always  present  in  digital  images  collected  using  digital  detectors  such  as  CCD  cameras. 
The  assumption  that  a  noise-free  image  is  the  convolution  of  the  true  object  and  the  PSF  means  that  the 
system  is  assumed  to  be  space  invariant.  For  many  applications,  this  is  a  valid  assumption;  however, 
even  when  the  system  violates  the  space- invariant  assumption,  this  model  still  produces  excellent  image 
restorations  as  will  be  demonstrated  in  Section  5.  We  call  the  PSF  invertible  if  its  Fourier  transform 
contains  no  zeros,  and  non-invertible  otherwise. 

The  object  and  PSF  estimates  are  generated  by  minimizing  a  cost  function  that  is  based  on  MLE 
theory  [34].  With  the  use  of  the  central  limit  theorem,  valid  for  all  but  quite  dim  objects,  we  model  nm(x) 
as  a  spatially-independent,  zero-mean  Gaussian  random  process  with  a  variance  equal  to  £t^(x)  +  /„,(x), 
where  o-^(x)  is  the  variance  of  the  Gaussian  noise  and  ;„,(x)  is  the  variance  of  the  Poisson  noise  in  the  m"' 
frame.  For  the  Gaussian  noise  model,  the  general  form  of  the  MLE-based  cost  function 

j[o(x),  A,  (x ),..., /j  j,  (x)]  that  is  minimized  to  generate  estimates  of  the  object  and  the  PSFs  is  given  by 


y  [(5(x),/j,  (x),...,  h„  (x)]  =  X  ,  [/,„  (x„ )-  /„  (x„  )f  +  yg[o(x„ ),  /i,  (x„ ),...,  h„  (x„ )],  (2) 

mini  V*n/ 
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where  M  is  the  number  of  measurement  frames,  N"  is  the  number  of  pixels  in  each  measurement  frame, 
o(x),/i,  are  the  estimated  object  and  PSFs,  /„{x)  is  the  estimated  m"’  measurement  frame,  g() 

is  a  function  used  for  regularization,  {x„}  are  discrete  spatial  locations  corresponding  to  pixels  in  a 
detector,  and  y  is  the  Lagrange  multiplier  weighting  the  relative  importance  of  the  regularization  term  and 
the  data-matching  term.  In  the  PCID  implementation  of  Eq.(2),  the  photon  noise  weighting  can  be  based 
on  the  measurement  as  shown,  or  can  be  based  on  the  estimate  of  the  measurement.  In  the  latter  case, 

/„,(x)  is  replaced  with  (x) . 

B.  PCID  implementation 

All  the  implementation  issues  for  the  PCID  algorithm  are  described  in  this  section.  These  issues  are:  (1) 
how  o(x)  and  /;„(x)  are  parameterized,  (2)  how  positivity  and  support  are  enforced,  (3)  how 

regularization  is  implemented,  (4)  the  means  that  we  use  to  minimize  j[6(x),/»|(x),...,/;„(x)]  ,  (5)  the 
method  to  achieve  de-aliasing-based  superresolution,  (6)  two  uniqueness  issues,  and  (7)  some  algorithm 
initialization  steps. 

1.  Parameterization  of  o(x)  and  A„(x) 

The  object  estimate  o(x)  is  always  parameterized  pixel  by  pixel  in  the  image  domain  so  that  we  do  not 
bias  the  solution  toward  preconceived  perceptions  of  what  it  should  be.  On  the  other  hand,  we  have 
implemented  two  choices  for  parameterizing  /i„,(x).  The  first  parameterization  is  by  its  image-domain 
pixel  values  as  for  the  object  estimate.  This  parameterization  is  the  most  general  one  but  can  make 
finding  global  minima  more  difficult  and  produces  restorations  with  higher  variances  than  achievable 
with  more  restrictive  parameterizations  (when  the  more  restrictive  parameterizations  are  valid).  The 
second  parameterization,  that  can  be  used  with  great  success  as  long  as  the  spectral  bandwidth  of  the 
measurements  isn’t  too  large,  is  in  terms  of  the  phase  values  of  the  field  in  the  pupil  of  the  imaging 
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system.  The  PSF  estimates  are  generated  from  the  pupil  phase  values  with  the  use  of  the  narrow-band 
Fourier  optics  model  [35].  We  choose  to  use  a  Zernike-based  modal  parameterization  [30],  but  others 
have  used  a  pi.\el-based  parameterization  [21]. 


2.  Positivity'  and  support 

Two  important  pieces  of  additional  knowledge  that  can  be  included  as  constraints  in  the  estimation 
process  are  the  finite  support  of  o(x),  if  applicable,  and  the  fact  that  the  intensities  of  o(x)  are  positive. 
The  support  constraint  is  implemented  simply  by  including  in  the  minimization  process  only  those  pixel 
locations  of  o(x)  that  are  in  the  known  support.  We  enforce  positivity  by  replacing  o(x)  in  Eq.(2)  with 
and  carrying  out  the  minimization  over  /?(x).  This  parameterization  of  o(x)  permits  the 
minimization  of  Eq.(2)  using  unconstrained  routines. 

3.  Regularization  approaches 

Two  versions  of  Eq.(2)  are  implemented  in  PCID  corresponding  to  the  two  methods  of  regularization. 
The  first  version  sets  y  =  0,  so  that  there  is  not  a  separate  regularization  term.  The  regularization  is 
implemented  in  the  creation  of  /„(x),  which  is  calculated  as 


^».(x)  =  /'mW*/».(x)*o(x),  (3) 

where  h^x)  is  a  regularization  filter.  A  benefit  of  carrying  out  regularization  as  shown  in  Eq.(3)  is  that 
well-understood  linear  filters  can  be  implemented.  The  potential  downside  of  this  manner  of 
regularization  is  that  the  problem  is  still  ill-conditioned  in  the  sense  that  the  global  minimum  of  Eq.(2)  is 
unaffected  by  spatial  frequency  values  of  o(x)  not  contained  in  the  Fourier-domain  support  of  /j,„(x).  We 
have  discovered  that  the  inclusion  of  hX\)  significantly  slows  the  amplification  of  these  spatial  frequency 
values  of  6(x),  so  that  regularization  is  easily  carried  out  by  restricting  the  number  of  iterations  the 
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minimization  procedure  is  allowed  to  execute.  We  note  that  the  final  estimate  of  o(x)  is  /7^(x)o(x),  not 
6(x).  A  second  benefit  of  this  type  of  regularization  that  follows  from  this  fact  is  that,  when  enforcing 
positivity  as  described  above,  the  resulting  unregularized  estimate  satisfies  the  positivity  constraint  but  the 
regularized  version  can  have  negatives  values  if  the  regularization  filter  has  negative  values.  Although 
the  positivity  constraint  can  reduce  noise  levels  in  images  [36]  and  improve  convergence  properties,  it 
also  limits  the  maximum  resolution  in  a  restored  object  to  less  than  is  possible  for  a  restoration  that  is 
permitted  to  have  negative  values  [37].  This  regularization  approach  maximizes  the  benefits  of  the 
positivity  constraint  while  removing  its  resolution  limitations.  The  image  restorations  produced  with  this 
approach  have  excellent  quality,  as  will  be  demonstrated  in  Section  5. 

The  second  version  of  Eq.(2)  that  is  implemented  in  PCID  uses  standard  Tikhonov  regularization 
[38];  thus,  g()  is  given  by 


g[o(x),/Ux),...,/J»  (x)]=  •  (4) 

n-1 

In  addition,  because  regularization  is  implemented  by  this  penalty  term,  the  estimated  measurement  is 
given  by 


<.W  =  *.(»)*'5(«)-  (5) 

A  key  benefit  of  using  Tikhonov  (or  other  penalty-term-based  regularization  schemes)  is  that  the  estimate 
o(x)  obtained  at  convergence  is  the  regularized  version,  so  it  is  unnecessary  to  regularize  by  iteration. 
This  property  is  useful  when  seeking  to  automate  the  algorithm  execution.  Our  primary  reason  for 
implementing  the  Tikhonov  penalty  term,  as  compared  to  other  penalty  terms,  is  to  make  calculating 
sample  variances  of  o(x)  to  compare  with  their  associated  CRBs  computationally  efficient. 
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Regardless  of  which  regularization  method  is  used,  the  regularization  parameter  must  be 
calculated  to  select  the  amount  of  regularization.  Many  approaches  to  this  calculation  have  been 
published  (see  Ref.  [39]  and  articles  referencing  it).  In  the  case  of  MFBD  applied  to  imaging  through 
atmospheric  turbulence,  a  straightforward  calculation  is  possible  as  is  now  described.  For  the  linear- 
filter-based  regularization  method,  we  choose  its  cutoff  frequency  (if  it  is  bandlimited)  or  its  half-width 
parameter  (if  it  is  not  bandlimited)  in  terms  of  where  the  signal-to-noise  ratio  (SNR)  of  the  azimuthally- 
averaged  image  energy  spectrum  estimate  is  one.  For  linear  image  restoration  using  speckle  imaging 
techniques,  it  has  been  shown  [40]  that  the  SNR  of  the  object  energy  spectrum  estimate  is  equal  to  the 
SNR  of  the  image  energy  spectrum  estimate.  Because  we  have  multiple  measurement  frames  in  MFBD, 
the  image  energy  spectrum  estimate  can  be  calculated  from  the  data  using  sample  statistics. 

For  the  second  regularization  approach,  the  image  energy  spectrum  SNR  is  also  calculated,  but 
choosing  and  implementing  the  regularization  parameter  is  more  challenging.  In  general.  Tikhonov 
regularization  produces  space-variant  regularization,  so  a  Fourier-domain  filter  representation  of  its 
regularization  is  not  valid.  Even  in  the  case  when  Tikhonov  regularization  is  space-invariant  (white 
Gaussian  noise),  knowledge  of  the  energy  spectra  of  the  blurring  PSFs  is  required.  We  currently  calculate 
the  regularization  parameter  to  provide  the  desired  ratio  of  the  data-matching  term  to  the  regularization 
term.  In  the  future,  we  plan  to  calculate  the  regularization  parameter  by  generating  an  analytical 
approximation  to  the  average  PSF  energy  spectrum  [40]  and  approximating  the  white  Gaussian  noise 
parameter  in  terms  of  the  measurement  noise  energy  spectrum  [9].  We  then  will  use  these  values  in  the 
white  Gaussian  noise  Tikhonov  filter  and  choose  its  regularization  parameter  to  produce  the  desired  half¬ 
width  dimension.  To  calculate  the  analytical  form  of  the  average  PSF  energy  spectrum,  the  Fried 
parameter  ro  must  be  estimated  from  the  blurred  imageiy'  (4 1  ]  or  from  ancillary  measurements. 

The  automatically-calculated  regularization  parameter  provides  a  good  starting  point  for  the 
image  restoration  process  and  is  especially  useful  for  algorithm  automation.  The  optimum  regularization 
parameter  is  a  function  of  how  the  restored  imagery  will  be  used;  for  this  reason,  we  include  the  option  of 
allowing  the  user  to  select  the  regularization  parameter  manually. 
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4.  Minimization  procedure 

We  use  a  conjugate  gradient  routine  to  minimize  the  cost  function  corresponding  to  the  regularization 
approach  that  is  implemented  and  all  the  additional  knowledge  that  is  applied.  We  provide  to  the 
minimization  routine  analytical  expressions  for  the  derivatives.  The  routine  that  we  use  is  based  on  the 
Numerical  Recipes  conjugate  gradient  routine  [42].  but  is  significantly  modified.  We  have  removed  the 
common  block  structure  to  make  the  routine  less  memory  intensive.  We  have  also  parallelized  the  routine 
as  will  be  discussed  in  Section  3.  We  have  implemented  four  modifications  to  help  speed  up  the 
minimization  process.  First,  we  choose  the  initial  step  sizes  and  directions  for  the  search  step  in  the 
routine  based  upon  the  MACOPT  approach  [43].  Second,  we  normalize  the  PSF  and  object  derivatives  to 
balance  the  impact  on  the  cost  of  mismatches  between  the  true  and  estimated  object  and  PSF  estimates. 
Third,  we  use  half-arrays  for  the  FFT  calculations  necessary  to  generate  i„,(x)  with  the  use  of  either 
Eq.(3)  or  Eq.(5).  Finally,  we  use  the  FFTW  software  package.  In  the  serial  version  of  PCID,  we  use  the 
2-D  version  of  FFTW,  but  in  the  parallelized  version  of  PCID  we  use  only  the  1-D  version  since  the  latest 
official  version  of  FFTW  (FFTW  3.1.2)  does  not  contain  a  distributed  memory  parallel  implementation. 
An  alpha  version  of  FFTW  3.2  is  available  that  does  have  a  distributed-memory  implementation,  and  we 
plan  to  evaluate  its  performance  compared  to  our  current  implementation  when  it  becomes  the  official 
version. 

5.  De-aliasing-based  superresolution 

It  is  not  unusual  to  have  measurement  data  that  is  aliased  due  to  size  of  the  object  being  imaged  relative  to 
the  number  of  pixels  in  the  observation  camera.  It  is  well  known  that  de-aliasing  can  be  accomplished  if 
multiple  sub-pixel-shifted  images  of  the  same  scene  are  available  [44];  however,  the  key  issue  in  the 
accuracy  of  the  de-aliased  Fourier  spectrum  is  the  accuracy  to  which  the  multiple  images  are  co-registered 
[45].  For  an  image  that  is  aliased  by  a  factor  of  X.  X‘  separate  and  linearly-independent  images  must  be 
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measured.  When  imaging  through  the  atmosphere  using  short-exposure  images,  hundreds  or  thousands  of 
images  can  be  collected,  so  this  requirement  is  easily  satisfied.  We  carr>’  out  the  de-aliasing  procedure  as 
part  of  the  MLE  estimation  and  jointly  estimate  the  measured  image  shifts  along  with  the  de-aliased 
object  and  PSF  restorations.  The  only  option  that  needs  to  be  selected  is  to  embed  the  measured  data 
arrays  in  larger  arrays  that  are  at  least  Nyquist  sampled. 

6.  Two  uniqueness  issues 

As  explained  in  [4],  there  are  two  requirements  for  the  support-constrained  blind  deconvolution  problem 
to  generate  unique  solutions.  The  first  is  that  the  sum  of  either  the  object  or  each  of  the  PSF  intensity 
values  must  be  held  constant.  The  reason  for  this  requirement  is  that  an  estimated  measurement  remains 
unchanged  if  the  object  estimate  is  multiplied  by  a  constant  and  the  PSF  estimate  is  divided  by  that  same 
constant.  For  pixel-based  blind  deconvolution,  it  is  necessary  then  to  formulate  the  estimate  of  the  m"* 
measurement,  for  all  m,  in  a  way  that  keeps  the  total  PSF  intensity  fixed.  We  accomplish  this  by 

replacing  /j„(x)  in  Eqs.  (3)  and  (5)  with  ^„(x)/^/»„(x).  This  forces  the  sum  of  the  object  intensity 

values  to  remain  fixed  and  thus  the  sum  of  the  intensity  values  for  each  PSF.  For  the  Zemike-based  PSF 
estimate  version  of  PCID,  the  sum  of  the  PSF  intensity  values  is  always  equal  to  one  because  of  the 
standard  mathematical  model  that  relates  the  Zemike  coefficients  and  polynomials  to  the  corresponding 
PSF  [35]. 

The  second  requirement  for  support-constrained  blind  deconvolution  to  be  unique  is  that  the 
locations  of  the  object  and  the  PSFs  must  be  fixed.  The  reason  for  this  requirement  is  that  an  estimated 
measurement  remains  unchanged  if  the  object  estimate  is  shified  a  certain  amount  in  a  given  direction  and 
the  PSF  estimate  is  shifted  the  same  amount  in  the  opposite  direction.  We  note  that  the  structures  of  both 
the  object  and  the  PSF  estimates  must  remain  unchanged  when  shifted  for  the  estimated  measurement  to 
remain  unchanged.  Since  the  PSFs  for  realistic  imaging  systems  have  infinite  support,  any  shift  in  a  PSF 
estimate  will  change  its  structure  because  of  wraparound  in  the  array  containing  the  PSF  estimate.  This 
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removes,  theoretically,  the  insensitivity  of  the  estimated  measurement  to  shifts  in  the  PSF  estimate.  For 
real  estimation  problems,  the  intensity  values  in  a  PSF  estimate  approach  zero  at  the  edges  of  the  array,  so 
the  blind  deconvolution  problem  is  ill-conditioned  if  the  object  and  PSF  locations  arc  not  fixed.  We 
remove  this  ill-conditioned  issue  by  applying  support  constraints  to  both  the  object  and  each  of  the  PSFs. 
To  minimize  the  deleterious  effects  of  applying  a  support  constraint  to  an  infinite-support  PSF,  we  choose 
the  size  of  the  support  to  encompass  at  least  99%  of  the  estimated  PSF  intensity 

7.  Initialization  procedures 

The  measurement  data  needs  to  be  preprocessed  to  remove  atmospheric  background,  camera  dark  current 
biases,  and  optical  system  flat  field  effects  prior  to  be  used  in  the  PCID  algorithm.  In  addition,  the 
measurement  data  must  be  scaled  so  that  the  intensity  values  correspond  to  photon  counts.  If  Zernike- 
based  PSF  estimation  is  used,  the  initial  estimates  of  the  Zemike  parameters  are  set  to  zero,  except  for  tip 
and  tilt  that  are  estimated  from  the  measurement  data.  For  pixel-based  PSF  estimation,  the  default 
initialization  is  an  Airy  pattern  corresponding  to  the  telescope  used  for  the  observations,  but  user- 
provided  initial  PSF  estimates  can  be  used  if  desired. 

An  initial  estimate  of  the  object  is  generated  by  co-adding  centered  versions  of  all  the 
measurement  frames.  The  object  support  region,  if  it  is  to  be  used  as  a  constraint,  is  generated  from  the 
initial  object  estimate  by  thresholding  the  initial  object  estimate  by  an  amount  calculated  from  the  camera 
read  noise  and  the  calculated  ro.  The  read  noise  is  used  to  make  sure  that  the  object  support  region  is  not 
extended  into  the  noise-only  parts  of  the  measurement  frame,  and  the  ro  parameter  is  used  to  estimate  the 
atmospheric  blurring  in  order  to  shrink  the  object  support  to  more  closely  match  the  true  object  support. 

An  important  feature  of  the  PCID  algorithm  is  that  it  can  use  clipped  frames  of  data  that  are 
brought  about  by  either  image  jitter  or  by  the  fact  that  the  object  is  too  large  for  the  imaging  system  field 
of  view.  As  long  as  all  parts  of  the  object  are  visible  when  considering  all  the  measurement  frames,  the 
PCID  algorithm  can  restore  a  single  undipped  estimate  since  it  includes  clipping  effects  in  the  forward 
model.  Clipping  is  a  problem,  however,  when  generating  either  an  initial  estimate  of  the  object  or  the 
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object  support.  For  this  reason,  we  have  implemented  an  automatic  way  to  detect  and  thus  to  exclude 
clipped  measurement  frames.  We  use  the  fact  that  the  energy  spectra  of  clipped  images  have  excessively- 
high  values  along  the  frequency  axis  corresponding  to  the  clipped  direction  to  automatically  identify 
clipped  images. 

3.  PCID  PARALLELIZATION 

We  describe  the  parallel  architecture  of  the  PCID  algorithm  in  this  section  and  discuss  the  features  of  the 
algorithm  that  ultimately  limit  its  scalability  across  multiple  processors.  We  then  demonstrate  its 
scalability  using  a  Cray  XDI  parallel  computer  that  is  representative  of  available  commodity  clusters, 
where  commodity  clusters  are  parallel  computers  that  use  non-special ized  hardware.  In  this  context,  the 
term  “scalability”  refers  to  the  ability  of  the  PCID  algorithm  to  use  increasing  numbers  of  processors 
efficiently  to  reduce  the  wall-clock  execution  time  [46].  For  example,  a  program  can  be  said  to  scale 
perfectly  if  the  wall-clock  execution  time  decreases  by  a  factor  of  two  when  the  number  of  processors 
used  increases  by  a  factor  of  two.  The  PCID  algorithm  is  executed  on  top  of  a  message-passing  interface 
(MPI)  environment  called  MPICH  [47]. 

A.  PCID  parallel  architecture 

The  PCID  algorithm  is  structured  as  a  multiple-instruction  multiple-data  (MIMD)  program  in  a  master- 
worker  architecture.  The  MIMD  structure  allows  commands  in  a  single  program  to  be  identified  as 
belonging  to  the  master  process  or  a  worker  process  and  only  executed  by  the  correct  process.  The  master 
process  controls  the  overall  program  execution  and  the  worker  processes  carry  out  all  the 
computationally-intensive  calculations.  A  top-level  flowchart  of  the  PCID  parallel  architecture  is  shown 
in  Fig.  I  where,  for  clarity,  only  one  worker  process  is  displayed.  We  next  describe  first  the  master 
process  and  then  the  worker  process. 

The  master  process  only  uses  one  processor  regardless  of  the  number  of  processors  available 
since  its  computational  load  is  much  lighter  than  for  the  worker  processes.  The  master  process  is 
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responsible  for  the  overall  algorithm  execution.  This  includes:  (1)  accepting  all  the  inputs  to  the  PCID 
algorithm,  (2)  reading  in  the  measurement  data  from  a  file,  (3)  calculating  all  the  algorithm  initialization 
values,  (4)  parceling  out  the  measurement  data,  the  processing  parameters,  and  the  initial  estimates  of  the 
object  and  the  PSFs  to  the  worker  processes,  (5)  receiving  back  from  the  worker  processes  the 
information  necessary  update  the  value  of  the  cost  function  and  to  select  the  next  conjugate  gradient 
direction,  (6)  passing  back  to  the  worker  processes  the  updates  to  the  estimates  of  the  object  and  PSFs, 
and  (7)  determining  when  to  halt  the  overall  program  execution  either  because  the  maximum  number  of 
iterations  has  been  reached  or  because  the  relative  change  in  the  cost  function  is  less  than  the  minimum 
value  selected.  When  Zemike-based  PSF  estimation  is  selected,  the  master  process  passes  to  each  worker 
process  the  information  that  it  needs  to  calculate  that  portion  of  the  Zemike  polynomials  it  will  use  in  the 
estimation  process. 

Each  worker  process  calculates  its  part  of  the  cost  function  and  the  gradient  of  the  cost  function 
with  respect  to  the  current  object  estimate  and  the  current  PSF  estimates  for  its  measurement  frames.  In 
addition,  each  worker  process  calculates  its  portion  of  the  quantities  needed  to  determine  the  next 
conjugate  direction.  The  calculations  carried  out  by  a  given  worker  process  are  independent  of  all  the 
other  worker  processes.  This  means  that  no  communications  are  necessary  between  worker  processes  and 
thus  the  only  communications  that  occur  are  between  the  master  process  and  each  of  the  worker 
processes.  When  each  worker  process  uses  only  one  processor,  the  communication  between  processors  is 
minimal  after  the  initialization  communications  and  thus  the  limiting  factors  in  PCID  execution  speed  are 
the  number  and  speed  of  the  processors.  For  this  reason,  the  benefit  of  this  level  of  parallelization  is  that 
it  scales  closely  with  the  number  of  processors  used  to  run  PCID.  As  a  result,  increasing  the  number  of 
processors  by  a  factor  of  K  decreases  the  execution  time  by  a  factor  of  1/K  (neglecting  the  initialization 
execution  time).  The  downside  of  this  level  of  parallelization  is  that  the  maximum  amount  of 
parallelization  possible  is  equal  to  the  number  of  measurement  frames.  For  further  increases  in  execution 
speed,  each  worker  process  must  use  multiple  processors.  Issues  associated  with  this  increased  level  of 
parallelization  are  discussed  next. 
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Fig.  1.  PCID  parallel  architecture 


Each  worker  process  is  required  to  carry  out  convolution  operations  in  the  course  of  calculating 
the  cost  function  and  cost  function  gradient.  Convolution  operations  are  carried  out  most  efficiently  in 
the  Fourier  domain.  For  this  reason,  the  2-D  FFT  plays  a  dominant  role  in  the  total  execution  time;  in 
fact,  profiling  results  show  that  over  half  of  the  execution  time  of  the  PCID  algorithm  is  spent  carrying 
out  2-D  FFTs.  For  this  reason,  it  is  essential  to  maximize  the  speed  of  the  2-D  FFT  operation.  The 
general  2-D  FFT  of  an  array  is  carried  out  by  first  replacing  its  rows  with  their  l-D  FFTs,  and  then 
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replacing  the  columns  of  this  array  with  their  l-D  FFTs.  When  only  one  processor  is  used  per  worker 
process,  all  the  data  needed  for  the  2-D  FFTs  are  contained  in  local  memory  and  thus  the  execution  speed 
of  the  2-D  FFT  is  driven  by  the  speed  of  serially  carry  ing  out  the  l-D  FFTs.  When  multiple  processors 
are  used  per  worker  process,  the  l-D  FFTs  can  be  carried  out  in  parallel  by  all  the  worker  processors. 
This  can  increase  the  execution  speed  of  the  2-D  FFT;  however,  data  transfer  is  necessary  between  the 
worker  processors  during  the  2-D  FFT  operation  which  mitigates  the  speed  improvements,  even  when 
high-speed  interconnects  such  as  Infiniband  [48J  and  Myrinet  [49]  are  employed.  For  this  reason,  we 
have  sought  to  minimize  the  communications  between  the  processors  in  a  worker  process  to  carry  out  the 
2-D  FFT.  Because  all  fast  l-D  FFTs  require  the  data  to  be  in  local  memory,  the  master  process  always 
spreads  the  data  across  worker  processors  by  row  so  that  the  first  set  of  l-D  FFTs  can  be  executed  with 
the  rows  currently  in  memory.  To  maximize  the  speed  of  execution,  each  worker  process  must  have  the 
same  number  of  rows.  After  the  l-D  row  FFTs  have  been  calculated,  the  rows  of  the  array  must  be 
redistributed  among  the  worker  processors  so  that  each  column  of  the  array  resides  in  the  local  memory  of 
a  worker  processor.  To  explain  how  we  do  this,  suppose  that  there  are  two  processors  in  a  worker  process 
as  shown  in  Fig.  2.  For  an  N  by  N  array  (N  even),  one  processor  is  sent  the  first  N/2  rows  and  the  second 
is  sent  the  last  N/2  rows.  After  l-D  FFTs  are  taken  of  each  row,  the  resulting  rows  are  split  into  two  N/2 
by  N/2  blocks  that  are  separately  transposed.  All  of  these  operations  take  place  in  local  memory,  so  no 
communication  is  necessary  between  the  worker  processors.  Next,  the  off-diagonal  blocks  are  swapped 
between  the  two  processors  -  this  is  the  only  step  that  requires  communication  between  the  two  worker 
processors.  Once  the  transfer  is  complete,  each  worker  processor  has  N/2  columns  of  the  array  and  can 
complete  the  2-D  FFT  in  local  memory.  For  the  general  case  of  W  worker  processors  per  worker  process, 
N  must  be  evenly  divisible  by  W.  Because  rows  and  columns  are  swapped,  communications  only  occur 
between  pairs  of  worker  processors,  permitting  the  use  of  point-to-point  MPI  communications. 
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Fig.  2.  PCID  parallel  2-D  FFT  flowchart 


Since  each  processor  must  have,  at  a  minimum,  an  entire  row  of  the  array  to  be  Fourier 
transformed  using  fast  2-D  FFT  routines,  the  maximum  number  of  processors  that  can  be  used  per 
measurement  frame  is  N.  When  M  measurement  frames  are  used  for  a  single  object  estimate,  the 
maximum  number  of  processors  that  PCID  can  use  is  MN.  One  example  of  how  many  processors  that 
can  be  used  by  PCID  is  for  the  case  of  arrays  that  are  1024  by  1024,  and  when  400  measurement  frames 
are  used,  as  is  the  situation  for  one  of  the  example  restorations  in  Section  5.  Theoretically,  409,600 
processors  can  be  used  in  the  restoration  process;  however,  intemode  communication  time  dominates  the 
algorithm  execution  time  for  far  fewer  processors. 
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B.  PCID  parallel  scalability  demonstration 

We  used  a  Cray  XDI  supercomputer  [33]  located  at  the  Maui  High-Performance  Computing  Center  [50] 
to  explore  the  scalability  of  PCID  as  a  function  of  the  number  of  worker  processors  used  in  a  worker 
process.  The  Cray  XDI  is  an  SLES  9.0  Linux-based  cluster  with  144  nodes,  each  with  two  2.4  GHz 
AMD64  Opteron  processors  and  8  GBytes  of  RAM.  Each  processor  has  a  theoretical  peak  speed  of -4.8 
GFLOPS.  The  nodes  are  connected  via  a  Fat  Tree  RapidArray  Interconnect  switch  that  has  a  bandwidth 
of  4  GBytes/second  with  a  MPI  latency  of  1.7  ps.  The  PCID  algorithm  uses  the  MPl  standard  to 
communicate  between  nodes.  For  maximum  speed,  it  uses  vendor-optimized  versions  on  the  platforms 
that  it  runs  on  when  available  (as  is  the  case  for  the  Cray  XDI).  but  is  fully  compatible  with  both  the 
MPICHI  and  MPICtl2  freely-available  implementations  of  MPI  [47].  Although  the  Cray  XDI  nodes 
have  two  processors  each,  we  chose  to  use  only  one  processor  per  node  for  the  scalability  tests  to  ensure 
that  all  interprocessor  communications  use  the  intemode  communication  channels,  not  on-board 
communication  channels.  This  is  the  severest  test  of  the  scalability  of  a  parallel  program. 

Execution  times  for  PCID  as  a  function  of  the  number  of  processors  included  in  a  worker  process 
for  measurement  frames  that  are  128  by  128  pixels  in  size  are  shown  in  Fig.  3a.  The  numbers  of 
processors  included  in  the  worker  process  were  limited  to  powers  of  two  because  the  array  size  is  a  power 
of  two.  These  results  were  generated  using  only  one  worker  process  and  letting  PCID  run  for  100 
iterations  (one  iteration  corresponds  to  a  single  conjugate-gradient  search  direction).  It  can  be  seen  that, 
as  the  number  of  processors  is  increased  from  1  to  32,  the  PCID  execution  times  decrease,  indicating  that 
the  execution  times  in  this  region  are  limited  by  the  availability  of  processor  cycles.  For  more  than  32 
processors,  the  execution  times  increase  as  the  number  of  processors  is  increased,  indicating  that 
internode  communications  are  limiting  the  execution  time.  To  quantify  the  scalability  of  PCID  across 
processors,  we  calculated  a  scalability  factor  for  each  number  of  processors  greater  than  one  using  the 
data  in  Fig.  3a.  For  a  given  number  of  processors,  say  2".  the  scalability  factor  is  equal  to  the  execution 
time  for  2"‘'  processors  divided  by  the  execution  time  for  2"  processors.  If  the  PCID  algorithm  scaled 
perfectly,  the  scalability  factor  would  be  equal  to  two  for  any  value  of  n.  If  adding  processors  did  not 
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help,  the  factor  would  be  one.  Values  of  the  factor  less  than  one  indicate  that  execution  times  are 
increasing  as  processors  are  added.  The  scalability  factor  is  plotted  in  Fig.  3b  as  a  function  of  the  number 
of  processors.  It  can  be  seen  that  PCID  only  scales  perfectly  when  moving  from  one  to  two  processors. 
Beyond  this,  the  scalability  factor  drops  to  around  1.5.  Finally,  for  more  than  32  processors,  the 
scalability  factor  is  less  than  one  because  execution  times  are  increasing. 


1  0  1  00 
Number  or  processors 


(a) 


Fig.  3.  (a)  Execution  times  and  (b)  scalability 

measurement  frames. 


(b) 


factors  for  the  PCID  algorithm  for  128  by  128 


We  next  generated  results  in  the  same  way  as  for  the  results  in  Fig.  3  except  we  used 
measurement  frames  that  are  512  by  512  pixels  in  size.  The  results  are  shown  in  Fig.  4.  It  can  be  seen 
that  the  scalability  for  larger  arrays  is  better  than  for  smaller  arrays.  This  is  due  to  the  fact  that,  for  larger 
arrays,  each  processor  has  more  data  in  local  memory  than  for  smaller  arrays,  thus  the  execution  time  is 
limited  by  the  CPU  cycles  available.  Eventually  the  intemode  communications  load  will  cause  execution 
times  to  increase  as  more  processors  are  added. 
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Fig.  4.  (a)  Execution  times  and  (b)  scalability  factors  for  the  PCID  algorithm  for  512  by  512 

measurement  frames. 


As  a  reminder,  the  results  shown  in  Figs.  3  and  4  demonstrate  the  scalability  of  adding  multiple 
processors  to  a  worker  process  for  a  single  measurement  frame.  If  there  are  multiple  measurement  frames 
to  be  included  in  the  restoration  process,  the  best  parallelization  approach  is  data  parallel  if  possible,  since 
the  scalability  is  essentially  perfect  in  this  case.  Only  when  there  are  more  processors  than  measurement 
frames  should  one  include  multiple  processors  in  a  worker  process. 

4.  CRB  THEORY  AND  RESULTS 

To  evaluate  the  quality  of  image  restorations  using  a  given  algorithm,  it  is  useful  to  have  algorithm- 
independent  limits  to  image  quality  as  benchmarks.  There  are  a  number  of  algorithm-independent  metrics 
that  can  be  used  to  evaluate  image  quality,  many  of  which  are  subjective  in  nature.  We  choose  to 
quantify  image  quality  in  terms  of  the  SNRs  of  the  image  restorations  in  order  to  have  an  objective  metric 
to  evaluate  algorithm  performance.  Cram6r-Rao  lower  bound  theoiy  provides  a  way  to  calculate 
algorithm-independent  lower  bounds  (CRBs)  to  the  variances  of  estimates  of  a  set  of  parameters  (thus 
upper  bounds  to  the  SNRs  of  parameter  estimates).  These  lower  bounds  can  be  calculated  for  unbiased  as 
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well  as  biased  estimates  of  the  parameters.  In  this  section  we  present  an  overview  of  CRB  theory  for 
unbiased  and  biased  estimates  of  a  set  of  parameters  and  then  apply  the  theory  to  unbiased  as  well  as 
regularized  (a  specific  type  of  bias)  estimates  of  the  image-domain  intensities  of  an  object.  In  particular, 
we  present  CRB  expressions  for  scenarios  that  we  also  can  calculate  meaningful  sample  variances  using 
PCID.  We  compare  these  CRBs  to  the  PCID  sample  variances  and  show  that  image  restorations 
produced  by  PCID  achieve  or  closely  approach  them. 

A,  CRB  theory 

Let  ®  be  a  vector  containing  the  parameters  to  be  estimated.  Any  unbiased  estimate  of  O,  d>„ ,  that  is 
based  on  a  set  of  measurements  y  has  variances  that  are  bounded  below  by  the  CRBs  that  are  the  diagonal 
elements  of  the  inverse  of  the  Fisher  information  matrix  (FIM)  F  corresponding  to  O  and  y. 
Mathematically,  we  have 

var(<i>„)>diag{F''),  (6) 

where  var(i„)  is  the  variance  of  and  diag(p')  denotes  the  diagonal  elements  of  F'.  The  element  of 
F  in  the  p*  row  and  q“'  column,  [Fjp,,  is  given  by 


[fL  =  e 


ain[/(y;a>)]ain|/(y;«l>)]1 

54»„  eo,  J’ 


(7) 


where  In  is  the  natural  logarithm,  f(y;«I>)  is  the  probability  density  function  of  the  measurement  y  that  is  a 
function  of  <D,  and  E{ }  is  the  expected  value  of  the  quantity  in  braces. 

Often,  biased  estimates  of  O  are  sought  because  either  unbiased  estimates  are  not  possible  or 
because  a  given  biased  estimate  has  more  desirable  properties  than  an  unbiased  estimate  such  as  lower 


mean-square  error.  Let  <i>^  be  a  biased  estimate  of  <b  with  bias  where  the  subscript  O  indicates  that 

the  bias  is  a  function  of  the  true  parameters.  Then  the  variances  of  <1>a  are  bounded  below  as  shown  in 
the  following  equation: 


var(a>J>diag[(l  +  V*b^)F-'(l+V^b^y]  .  (8) 

where  V^b^  is  the  gradient  of  b«  with  respect  to  (I>  and  I  is  an  identity  matrix  of  the  appropriate 
dimension. 

To  apply  CRB  theory  to  the  MFBD  image  restoration  problem,  the  imaging  model  of  Eq.(l)  must 
be  rewritten  in  a  discrete  vector  form.  To  this  end,  let  a  be  a  vector  that  contains  the  spatial  locations  of 
the  intensity  values  of  i„ix)  as  would  be  obtained  for  an  image  collected  with  a  CCD  array.  Then  let  y^, 
0,  and  q„,  be  one-dimensional  vectors  that  contain  the  values  of  /„,(«).  o(ct),  and  respectively,  on 
the  grid  defined  by  a.  In  addition,  let  Hm  be  the  block-circulant  system  matrix  associated  with  h„,(a)  [5 1  ]. 
This  permits  rewriting  Eq.(l)  as  a  matri.x-vector  equation  given  by 

y„,  =  H„o+n„.  (9) 

Now  let  all  M  measurements  be  concatenated  into  a  single  measurement  vector  y,  where  y  =  [y . ,yMr- 

Then,  using  Eq.(9),  we  can  write  the  concatenated  measurement  equation  as 

y  =  He  +  n,  (10) 
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where  H0  =  [H|0 . Hm0]^.  t)  =  [r|i,...,  t|m]^  and  the  PDF  of  t|  is  denoted  by  f,,(r|).  Because  both  photon 

and  camera  read  noises  are  statistically  independent  across  pixels  and  across  measurement  frames,  f,,(Ti)  is 
equal  to  the  multiplication  of  all  the  single-pixel  PDFs. 

Using  Eqs.  (6)  -  (8)  and  (10),  we  can  derive  CRB  expressions  for  unbiased  and  biased  estimates 
of  restored  images  that  can  be  compared  to  the  corresponding  PCID  restored-image  sample  variances; 
however,  prior  to  doing  so,  two  issues  must  be  addressed.  The  first  issue  is  that  an  unbiased  estimate  of 
the  underlying  true  object  must  be  attainable.  For  a  realistic  imaging  system,  an  unbiased  estimate  is  not 
attainable  because  its  PSF  is  not  invertible.  Theoretically,  an  unbiased  estimate  is  attainable  when 
utilizing  a  support  constraint  as  a  result  of  superresolution  even  with  a  non-invertible  PSF;  however,  with 
any  noise  at  all,  all  superresolved  Fourier  spectra  are  biased  [52].  Nonetheless,  we  will  derive  unbiased 
image-restoration  CRBs  for  two  reasons.  The  first  is  that  we  can  employ  invertible  PSFs  in  order  to  test 
the  performance  of  the  PCID  algorithm.  The  second  reason  is  that  we  can  use  a  particular  form  of  the 
FIM  for  non-blind  unbiased  CRBs  to  calculate  non-blind  regularized  CRBs  for  non-invertible  PSFs  for 
regularizing  filters  that  are  admissible  (as  defined  in  the  next  paragraph). 

The  second  issue  is  the  calculation  of  biased  CRBs.  It  can  be  seen  from  Eq.  (8)  that  an  analytical 
expression  for  the  gradient  of  the  bias  must  be  available  to  calculate  biased  CRBs.  Such  an  expression  is 
available  for  regularized  image  restorations  when  the  Fourier-domain  support  of  the  regularization  filter  is 
contained  in  the  Fourier-domain  support  of  the  Fourier-transformed  PSF  but  is  generally  not  available 
otherwise.  (We  will  refer  to  this  type  of  regularization  filter  as  an  admissible  regularization  filter).  It  is 
not  available  for  any  type  of  regularization  that  includes  estimated  Fourier  spectra  outside  the  passband  of 
the  PSF,  including  all  penalty-based  regularization  methods  such  as  Tikhonov  regularization.  The  reason 
that  analytical  expressions  are  not  available  for  non-admissible  regularization  is  because  this  type  of 
regularization  passes  restored  frequency  content  outside  the  bandwidth  of  the  measurement  system;  that 
is,  they  pass  superresolved  frequency  content.  As  stated  in  the  previous  paragraph,  all  superresolved 
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Fourier  data  are  biased,  and  analytical  expressions  for  the  biases  are  not  known  except  in  limited 
situations. 

Another  situation  where  bias  gradient  expressions  are  not  available  is  when  inequality  constraints 
such  as  positivity  are  used  in  the  image  restoration  process.  Initial  steps  have  been  taken  to  generate  an 
alternate  way  to  calculate  inequality-constrained  biased  CRBs  [36],  but  much  more  work  needs  to  be 
carried  out  before  this  approach  can  be  applied  to  positivity-constrained  image  restoration. 

Because  of  these  two  issues,  there  are  only  a  limited  number  of  scenarios  for  which  both  PCID 
image  restorations  can  be  generated  and  CRBs  can  be  calculated.  A  support  constraint  can  always  be 
included  in  any  of  these  scenarios.  All  scenarios  that  include  unbiased  image  restorations  fit  in  this 
category.  Almost  no  scenarios  that  include  biased  image  restorations  fit  in  this  category  because  of  the 
two  ways  that  PCID  implements  regularization.  The  Fourier-filter-based  regularization  approach 
accomplishes  the  requisite  regularization  by  limiting  the  number  of  iterations.  This  produces  biases  in  the 
image  restorations  that  can  not  be  calculated  analytically.  The  Tikhonov  regularization  approach  allows 
the  algorithm  to  iterate  to  convergence;  however,  the  image  restorations  contain  spatial  frequency  content 
from  the  entire  Fourier  domain.  As  a  result,  only  scenarios  that  include  invertible  PSFs  have  biases  that 
can  be  calculated  analytically.  In  addition,  our  CRB  calculations  assume  space-invariant  regularization, 
so  only  white  measurement  noise  can  be  modeled  when  Tikhonov  regularization  is  employed  [38]. 

As  a  result  of  the  issues  discussed  above,  we  have  been  able  to  compare  sample  variances  of 
PCID  image  restorations  to  their  associated  CRBS  for  the  following  scenarios:  (I)  non-blind  unbiased 
image  restorations  for  any  combination  of  read  and  photon  noises,  (2)  non-blind  biased  image  restorations 
for  invertible  PSFs,  Tikhonov  regularization,  and  white  read  noise,  and  (3)  blind  unbiased  image 
restorations  for  pixel-based  invertible  PSF  estimation  for  any  combination  of  read  and  photon  noises.  No 
blind  image  restorations  that  use  Zemike-based  PSF  estimation  are  included  in  the  scenarios  because  all 
Zemike-based  PSFs  are  non-invertible.  Any  support  constraints  that  contain  the  true  support  regions  for 
the  object  and/or  PSFs  can  be  included  in  any  of  these  scenarios. 
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Next,  we  present  an  expression  for  the  FIM  to  be  used  in  Eqs.  (6)  and  (8)  to  generate  the  CRBs 
that  will  be  compared  to  PCID  sample  variances.  We  present  the  FIM  for  blind  deconvolution  with  pixel- 
based  PSF  estimation  since  the  FIM  for  non-blind  deconvolution  is  a  special  case  of  the  blind  FIM.  In 
blind  deconvolution,  the  functionally-independent  parameters  that  generate  the  elements  of  both  {H,„} 
and  0  are  estimated,  where  {H^}  denotes  the  set  H„„  m=l,...,M.  Let  O  be  the  vector  that  contains  all  of 

N 

these  parameters.  For  pixel-based  estimation  of  /i„,(a)  and  o(a)  when  ^ /?,„(«„)  is  kept  fixed,  only  N-1 

n=I 

elements  of  each  of  the  h„{a)  are  functionally  independent.  Without  loss  of  generality,  we  assume  that 
the  first  N-1  elements  of  each  of  the  /?m(a)  are  used  to  generate  each  corresponding  H,„  and  thus  are  the 
elements  of  the  vector  C*  along  with  the  N  elements  of  0.  The  FIM,  F,  associated  with  y,  O,  and  q  has  a 
block  structure  that  is  given  by 


F  = 


^11  •  .  • 


(11) 


where  F„,i  =  Fim  (i.e.,  the  FIM  is  symmetric).  In  terms  of  this  block  structure,  the  first  column,  the  first 
row,  and  the  diagonal  are  the  only  parts  of  F  that  are  non-zero.  Using  the  results  in  [53],  it  can  be  shown 
that  the  element  of  Fn  in  the  p"'  row  and  the  q""  column,  [Fii]pq,  is  given  by 


tr.  1 

1*^11  Jpq  “ZjZj  .  /  \  2(  \ 


m=l  n=l 


(12) 


where  <t„^,(x)  is  the  variance  of  the  read  noise  in  the  m"'  image.  Similarly,  the  element  of  F,,,!  in  the  p"* 
row  and  the  q"'  column,  [Fn,i]pq,  is  given  by 
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(13) 


Fp  1  ^  K  (x„  -  X;.  )[»(x,.  -  X,  )-  o(x„  -  X  ,.>  )] 


Finally,  the  element  of  F,„m  in  the  p“’  row  and  the  q"'  column,  [F„,m]pq,  m  £  2,  is  given  by 


rp  1  ^  [o(x,.  -xJ-o(x„  -x  ^,)Iq(x„  -xJ-o(x„  -x^.)] 


(14) 


Equations  (12)  -  (14)  give  the  elements  of  F  without  support  constraints.  To  include  support  constraints, 
the  elements  of  O  are  restricted  to  contain  only  the  locations  of  {h„,{a)}  and  o(a)  that  are  inside  the 
support  constraint  regions,  where  {/»»,(«)}  denotes  the  set  /?„,(a),  for  m=l,...,M.  In  addition,  the  non¬ 
blind  version  of  F  is  equal  to  [Fn]. 

The  above  expression  for  F  can  be  used  to  generate  the  CRBs  for  unbiased  estimates  of  the  object 
and  PSF  intensities  with  the  use  of  Eq.  (6).  To  generate  regularized  CRBs,  the  bias  gradient  term, 
I  +  Vjibj,,  must  be  calculated  for  the  regularizing  filter  that  will  be  used.  The  general  form  of  the  bias 
gradient  term  is  block  diagonal  as  shown  below: 


I  +'^®bq, 


(15) 


where  the  Gn  block  is  used  to  generate  the  regularized  image  CRBs,  and  the  Gmm  block,  m  >  2,  is  used  to 
generate  regularized  CRBs  for  the  (m-l)'*'  PSF.  Our  interest  is  in  regularized  image  restorations,  not 
regularized  PSF  restorations.  For  this  reason,  we  set  all  blocks  of  the  bias  gradient  term  to  zero  except  for 
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Gii.  For  a  specified  regularizing  PSF  the  element  of  Gn  in  the  p'*'  row  and  the  q"'  column.  [GnJpq. 

is  given  by 


[GuL  (16) 

B.  Comparison  ofPCID  sample  variances  and  CRBs 

In  this  section  we  present  the  results  of  comparing  sample  variances  of  parameter  estimates  generated 
using  PCID  to  their  associated  CRBs  to  show  how  closely  the  PCID  sample  variances  approach  their 
theoretical  limits.  For  unbiased  image  restorations,  the  parameters  are  the  true  object  intensity  values  for 
each  pixel  in  the  object  support.  For  regularized  image  restorations,  the  parameters  are  the  regularized 
object  intensity  values  in  the  object  support.  Other  sets  of  parameters  could  have  been  chosen;  for 
example,  the  comparison  could  have  been  made  in  terms  of  the  power  spectra  of  the  restored  imagery. 
We  chose  to  carry  out  the  comparison  in  terms  of  intensity  values  in  the  restored  imagery  because  this  is  a 
useful  and  convenient  metric.  Since  functional  transformations  of  parameters  have  corresponding 
functional  transformations  of  their  CRBs  [25],  the  selection  of  a  specific  parameter  set  will  change  the 
absolute  values  of  the  sample  variances  and  the  CRBs,  but  not  the  results  of  the  comparison.  The 
comparisons  in  this  section  are  carried  out  for  the  three  scenarios  described  in  Section  4.A.:  ( 1 )  non-blind 
unbiased  image  restorations  corrupted  by  photon  and/or  read  noises,  (2)  non-blind  biased  image 
restorations  for  invertible  PSFs,  Tikhonov  regularization,  and  white  read  noise,  and  (3)  blind  unbiased 
image  restorations  corrupted  by  photon  and/or  read  noises  for  pixel-based  invertible  PSF  estimation.  We 
show  that  the  PCID  image  restorations  have  sample  variances  that  either  achieve  or  closely  approach  their 
associated  CRBs. 

To  carry  out  the  sample  variance/CRB  comparison  for  a  given  imaging  and  image  restoration 
scenario,  we  first  generated  50  measurement  data  sets  that  differed  only  in  their  realizations  of  the  read 
and  photon  noise  values.  Most  of  the  measurement  data  sets  consisted  of  one  frame  of  data,  but  a  few  of 


them  contained  ten  frames  of  data.  We  then  used  PCID  to  generate  50  restored  images,  one  for  each 
measurement  data  set..  Next  we  calculated  sample  variances  and  biases,  pixel  by  pixel,  using  the  50 
restorations.  For  unbiased  restorations,  we  calculated  the  biases  in  terms  of  the  true  object  being  imaged. 
For  regularization-based  biased  restorations,  we  calculated  the  biases  in  terms  of  the  regularized  version 
of  the  true  object.  Because  of  the  way  that  we  calculated  the  biases,  they  should  be  negligible  relative  to 
the  variances.  We  used  this  fact  to  ensure  that  the  PCID  image  restorations  correspond  to  the  true  global- 
minimum  solutions.  Once  this  was  verified,  we  calculated  the  ratio  of  the  sum  of  the  sample  variances  to 
the  sum  of  the  CRBs.  We  refer  to  this  ratio  as  the  PCID/CRB  ratio.  Numbers  close  to  one  show  that  the 
PCID  image  restorations  are  achieving  the  CRBs.  Because  of  the  large  number  of  comparisons  that 
underlie  the  results  in  this  section,  only  the  minimum,  maximum,  and  mean  values  of  this  ratio  are 
presented  for  each  of  the  three  scenarios  described  above. 

The  two  true  objects  used  in  the  comparison  study,  shown  in  Fig.  5,  are  a  computer-simulated 
satellite  model  called  OCNR  and  a  two-disk  object  called  two-circ.  The  satellite  object  was  chosen 
because  of  its  high  dynamic  range  and  significant  level  of  detail.  The  two-disk  object  was  chosen  to  be  a 
counterpoint  to  the  satellite  model  by  having  only  two  intensity  levels  (one  in  each  disk)  and  a  simple 
structure. 


Fig.  5.  The  two  objects  used  in  the  PCID  sample  variance  and  CRB  comparison:  (a)  OCNR,  (b)  two- 


Two  PSFs  were  used  for  the  comparison,  both  invertible  as  required  for  the  comparison.  The  first 
PSF,  called  tripsf,  is  the  inverse  Fourier  transform  of  the  following  function  H(f): 

f/(f)=  (1  -0.97f,  //.  Xl  -0.97f,  //J,  (17) 

where  f  =  (fi,f2)  and  f^  is  the  maximum  spatial  frequency  in  one  dimension  in  the  arrays  containing  the 
true  objects.  Because  H(f)  has  no  zeros,  it  produces  an  invertible  PSF.  The  second  PSF,  called  atmpsf, 
was  created  as  follows.  First,  we  generated  a  computer-simulated  image  of  an  unresolved  star  as  seen 
through  turbulence  and  an  imaging  system  characterized  by  the  ratio  D/to  =  8,  where  D  is  the  diameter  of 
the  imaging  system  [9].  This  image  is  a  PSF,  but  it  is  non-invertible  because  the  imaging  system  was 
modeled  realistically.  We  converted  it  to  being  invertible  by  multiplying  it  by  a  circular  support  region 
with  a  radius  of  10  pixels  that  contained  99%  of  the  true  PSF’s  energy.  For  some  of  the  image 
restorations,  we  used  ten  measurement  frames  with  this  style  of  turbulence  PSF.  In  this  situation,  ten 
separate  PSFs  were  generated,  all  with  the  same  D/r„  ratio,  but  with  separate  atmospheric  turbulence 
realizations.  The  PSF  corresponding  to  Eq.  (17)  and  one  of  the  turbulence  PSFs  are  shown  in  Fig.  6. 

The  three  object  support  constraint  options  available  for  use  by  PCID  were  the  true  object 
support,  a  larger  version  of  the  true  object  support  created  by  blurring  the  true  support  with  a  2  by  2 
blurring  kernel  (blur2  support),  and  the  smallest  centered  circular  support  region  that  contained  the  true 
object  (circle  support).  When  canying  out  blind  deconvolution  using  PCID,  only  the  atmpsf  PSFs  were 
used.  For  these  PSFs,  the  three  PSF  support  constraint  options  were  the  true  support,  a  centered  circle  of 
radius  20  pixels,  and  a  centered  circle  of  radius  30  pixels. 

For  scenario  I,  the  two  noise  options  were  read  noise  only  and  photon  plus  read  noise.  The  latter 
noise  model  was  dominated  by  the  photon  noise  component  in  total,  but  the  read  noise  component  was 
included  to  eliminate  the  possibility  of  dividing  by  zero  in  pixels  where  the  photon  noise  component  was 
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zero  (see  Eq.  (12)).  For  scenario  2,  only  read  noise  was  modeled  for  the  reasons  given  previously.  For 
scenario  3,  both  noise  options  were  used. 

The  regularized  results  in  scenario  (2)  used  a  value  of  the  Tikhonov  weighting  parameter  y  =  0. 1 . 
The  Fourier-domain  filter  corresponding  to  white  Gaussian  noise  with  the  measurement  PSF  and  this 
weighting  parameter  was  generated  using  the  expression  in  [38].  The  use  of  this  expression  required  that 
only  one  measurement  frame  be  included  in  each  image  restoration. 


m 


(a)  (b) 

Fig.  6.  The  two  PSFs  used  in  the  PCID  sample  variance  and  CRB  comparison  study:  (a)  tripsf,  (b)  one 
realization  of  atmpsf.  These  images  are  magnified  in  size  by  a  factor  of  8  as  compared  to  the  images 
shown  in  Fig.  5  in  order  to  better  display  their  detail. 


The  minimum,  maximum,  and  mean  values  of  the  PCID/CRB  ratios  for  all  three  scenarios  are 
shown  in  Table  1.  Theoretically,  the  values  of  the  PCID/CRB  ratios  should  be  no  less  than  one  for  any 
image  restoration;  however,  due  to  statistical  variations  in  the  sample  variances,  some  of  the  ratios  are 
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slightly  less  than  one.  The  scenario  1  results  were  generated  from  24  separate  combinations  of  objects, 
object  supports,  PSFs,  and  types  of  noise.  All  the  results  used  one  measurement  frame.  The  PCID/CRB 
ratios  most  closely  clustered  around  one  for  scenario  1,  clearly  showing  that  PCID  achieved  the  CRBs  for 
all  the  image  restorations.  The  scenario  2  results  were  generated  using  the  read-noise-only  subset  of  12 
image  restorations  used  for  the  scenario  1  results.  The  spread  of  the  PCID/CRB  ratios  is  greater  than  for 
scenario  1,  but  the  mean  is  still  essentially  one.  The  scenario  3  results  were  generated  using  only  the 
atmpsfs,  ten  measurement  frames  per  image  restoration,  and  the  true  and  blur2  object  supports.  Because 
scenario  3  is  the  blind-deconvolution  scenario,  the  PSF  support  constraints  described  above  were  also 
used.  Because  either  the  object  support  constraint  or  the  PSF  support  constraint  must  be  the  true  support 
to  keep  the  blind  deconvolution  problem  unique,  image  restorations  had  at  most  one  support  constraint 
larger  than  the  true  support.  Sixteen  combinations  of  these  options  were  used  to  generate  the  scenario  3 
results  in  Table  1.  The  PCID  sample  variances  are  higher  than  for  scenarios  1  and  2,  but  they  still  are 
only  8%  higher  than  the  CRBs 

Table  1 .  Mean,  minimum,  and  maximum  values  of  the  PCID/CRB  ratios  for  the  scenarios  described  in 

the  text. 


Scenario 

PCID/CRB  Ratios 

Mean 

Minimum 

Maximum 

1 

0.9854 

0.9402 

1.0326 

2 

1.0135 

0.9499 

1.0732 

3 

1 .0804 

0.9791 

1.205 

To  achieve  the  CRBs,  the  global  minimum  of  the  cost  function  must  be  found.  As  the  noise 
levels  increase,  and  as  the  atmospheric  blurring  increases,  it  becomes  more  difficult  to  find  the  global 


minimum  and  thus  to  achieve  the  CRBs.  We  are  working  on  ways  to  increase  the  robustness  of  PCID  in 
these  stressing  scenarios. 

5.  IMAGE  RESTORATION  EXAMPLES 

In  this  section  we  show  image  restorations  of  two  space  objects  to  demonstrate  the  quality  of  restorations 
possible  with  PCID.  Both  objects  were  imaged  using  telescopes  that  are  part  of  the  Maui  Space 
Surveillance  System  (MSSS)  [54].  The  two  objects,  the  Space  Shuttle  and  the  Hubble  Space  Telescope 
(HST),  were  chosen  for  the  image  restorations  because  they  are  familiar  to  a  wide  audience  and  because 
good  detail  can  be  seen  in  the  restorations.  The  use  of  multiple  images  of  the  same  object  increases  both 
the  restored  image  quality  and  the  probability  that  the  PCID  algorithm  will  find  the  global  minimum  of 
the  cost  function  shown  in  Eq.  (2);  however,  the  unblurred  object  has  to  be  the  same  in  all  the 
measurement  frames.  For  space  objects  orbiting  the  earth  at  low  altitudes  such  as  the  Space  Shuttle  and 
the  HST,  typical  observation  times  for  which  the  aspect  ratio  of  the  object  remains  unchanged  is  just  a 
few  seconds.  The  restorations  presented  in  this  section  use  all  the  measurement  frames  available  within 
the  allowed  observation  times. 

S.I  Space  Shuttle  image  restorations 

The  measurement  data  for  the  Space  Shuttle  restorations  were  collected  with  a  camera  mounted  on  the 
MSSS  1 .6  meter  telescope  without  the  use  of  adaptive  optics  compensation.  Representative  measurement 
frames  are  shown  in  Fig.  7.  The  camera  has  a  silicon-based  CCD  detector  that  is  128  by  128  pixels  in 
size,  four  readout  channels  to  increase  the  frame  rate,  and  a  read  noise  of  4  electrons.  The  spectral  filter 
through  which  the  light  was  collected  is  120  nm  in  width  with  a  mean  wavelength  of  940  nm.  The 
exposure  time  for  each  measurement  frame  was  0.9  msec,  ensuring  that  the  atmospheric  turbulence 
blurring  was  essentially  fixed  in  each  frame.  The  camera  frame  rate  was  250  Hz,  so  several  hundred 
measurement  frames  were  able  to  be  included  in  a  single  image  restoration.  Because  of  the  large  angular 
extent  of  the  Space  Shuttle,  the  telescope’s  field  of  view  (FOV)  was  set  to  1 20  microradians.  At  the  mean 
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observation  wavelength  of  940  nm,  the  camera’s  sampling  rate  was  only  30%  of  the  Nyquist  rate, 
resulting  in  data  that  was  aliased  by  a  factor  of  almost  four.  In  addition,  the  Space  Shuttle  was  slightly 
larger  than  this  FOV  even  when  perfectly  centered.  Due  to  residual  atmosphere  and  telescope  jitter,  many 
of  the  raw  images  are  clipped  a  significant  amount.  It  will  be  described  later  in  this  section  how  these 
negative  effects  are  accounted  for  in  the  image  restoration  process. 

Another  issue  associated  with  the  angular  size  of  the  Space  Shuttle  is  that  it  is  many  times  larger 
than  the  isoplanatic  angle  (the  angular  extent  over  which  atmospheric  turbulence  can  be  modeled  as  space 
invariant).  The  PCID  forward  model  assumes  that  the  image  can  be  modeled  as  a  convolution  of  the 
atmospheric  blurring  and  the  unblurred  object,  an  assumption  not  satisfied  by  the  data.  As  will  be  shown, 
though,  the  image  restorations  are  excellent  and  contain  near-diffraction-limited  detail.  This  indicates 
that  the  PCID  algorithm  is  robust  with  respect  to  violations  of  the  isoplanatic  assumption  governing  the 
forward  model  of  the  measurement. 

To  account  for  the  clipping  artifacts  in  the  measurement  frames,  the  128  by  128  measurement 
frame  arrays  were  embedded  in  the  center  of  256  by  256  arrays  with  zero  values  in  the  extra  pixels.  The 
embedded  measurement  frames  now  have  a  FOV  of  240  microradians.  The  embedding  process,  of 
course,  does  not  remove  the  clipping  in  the  measurements.  What  it  does  do  is  permit  the  restoration  of  an 
undipped  version  of  the  Space  Shuttle  in  a  256  by  256  array  since  the  forward  model  in  PCID  relating  the 
object  estimate  and  each  PSF  estimate  to  the  measured  data  incorporates  the  1 20  microradian  FOV  and  the 
size  of  the  Space  Shuttle  easily  fits  in  a  256  by  256  array. 

To  make  possible  unaliased  image  restorations,  each  embedded  measurement  frame  was 
upsampled  from  256  by  256  pixels  to  1024  by  1024  pixels  by  representing  a  single  pixel  in  the  256  by 
256  array  with  a  4  by  4  block  in  the  1 024  by  1 024  array,  where  each  value  in  the  4  by  4  block  is  the  value 
in  the  single  pixel.  Each  upsampled  measurement  frame  was  then  spatially  smoothed  with  a  4  by  4 
smoothing  kernel.  For  this  amount  of  aliasing,  a  minimum  of  16  iinearly-independcnt  data  frames  are 
needed  to  accomplish  the  de-aliasing,  a  limit  exceeded  by  more  than  a  factor  of  ten  for  the  restorations  in 
this  section. 
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Fig.  7.  Representative  measurement  frames  of  the  Space  Shuttle. 


(a) 


(b) 


Fig.  8.  Two  image  restorations  of  the  Space  Shuttle  produced  by  the  PCID  algorithm. 


The  PCID  algorithm  was  configured  to  estimate  the  PSFs  in  each  measurement  frame  using  the 
Zern ike-polynomial  parameterization  option.  The  first  one  hundred  Zemike  polynomials  were  used  in  the 
estimation  process.  All  of  the  Zernikc  coefficients  were  initialized  to  zero  except  for  the  tip  and  tilt 
values  that  were  initialized  for  each  PSF  to  the  tip  and  tilt  values  in  the  corresponding  measurement 
frame.  A  support  constraint  for  the  object  restoration  was  generated  automatically  as  described  in  Section 
2  and  used  in  the  restoration  process  as  was  a  positivity  constraint. 

Two  restored  images  of  the  Space  Shuttle  from  two  separate  parts  of  the  pass  are  shown  in  Fig.  8. 
The  restoration  in  Fig.  8a  was  achieved  with  the  use  of  400  measurement  frames.  Notice  the  high  quality 
of  the  restoration  as  demonstrated  by  the  detail  seen  in  the  cargo  bay  and  on  the  nose.  In  particular,  the 
two  round  windows  right  next  to  the  cargo  bay  are  resolved  to  near  the  diffraction  limit  of  the  optical 
system.  The  restored  image  shows  all  parts  of  the  Space  Shuttle  with  no  clipping  because  every  part  of 
the  Space  Shuttle  is  visible  in  at  least  one  measurement  frame  used  in  the  restoration  process.  The 
restored  image  of  the  shuttle  in  Fig.  8b  was  created  using  300  measurement  frames.  It  is  of  lesser  quality 
because  of  the  fewer  number  of  frames.  Notice  that  the  tip  of  the  nose  is  clipped  in  the  restoration  as  a 
result  of  it  being  clipped  in  all  the  measurement  frames  used  for  the  restoration.  Even  though  the  nose  is 
clipped,  the  rest  of  the  restored  image  is  unaffected  by  the  clipping  artifact. 

5.2  HST  image  restorations 

The  measurement  data  for  the  HST  restorations  were  collected  using  the  MSSS  A  EOS  3.6  meter 
telescope  with  its  visible-spectrum  adaptive  optics  system  [55].  Two  representative  measurement  data 
frames  are  shown  in  Fig.  9.  The  adaptive  optics  system  has  a  941  actuator  deformable  mirror  and  a 
Shack-Hartmann-based  wavefront  sensor.  The  camera  used  to  record  the  data  has  a  silicon-based  CCD 
detector  that  is  5 1 2  by  5 12  pixels  in  size,  only  one  readout  channel,  and  a  read  noise  of  1 2  electrons.  The 
spectral  filter  through  which  the  light  was  collected  is  360  nm  in  width  with  a  mean  wavelength  of  880 
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nm.  The  exposure  time  for  each  measurement  frame  was  98  msec.  The  exposure  time  was  much  greater 
than  the  atmospheric  correlation  time;  however,  because  of  the  adaptive  optics  correction,  the  effective 
system  correlation  time  was  much  greater  than  the  atmospheric  correlation  time.  The  camera  frame  rate 
was  2.4  Hz,  so  that  only  six  measurement  frames  were  included  in  a  single  image  restoration.  The  FOV 
of  the  telescope  was  set  to  5 1  microradians,  large  enough  to  contain  the  HST  and  small  enough  to  produce 
Nyquist-sampled  data.  In  addition,  the  degree  of  anisoplanicity  was  much  less  than  for  the  Space  Shuttle 
image  restorations. 

As  for  the  Space  Shuttle  image  restorations,  the  PCID  algorithm  was  configured  to  estimate  the 
PSFs  in  each  image  using  the  Zemike-polynomial  parameterization  option.  The  first  one  hundred  Zernike 
polynomials  were  used  in  the  estimation  process.  All  of  the  Zernike  coefficients  were  initialized  to  zero 
except  for  the  tip  and  tilt  values  that  were  initialized  for  each  PSF  to  the  tip  and  tilt  values  in  the 
corresponding  measurement  frame.  A  support  constraint  for  the  object  restoration  was  generated 
automatically  as  described  in  Section  2  and  used  in  the  restoration  process  as  was  a  positivity  constraint. 


Fig.  9.  Two  representative  measurement  frames  of  the  HST. 


(a)  (b) 


Fig.  10.  Two  PCID  restorations  of  the  HST. 


Two  HST  restorations  are  shown  in  Fig.  10.  There  is  a  noticeable  increase  in  image  quality  in 
these  restorations  when  compared  to  the  representative  measurement  frames;  however,  the  amount  of 
improvement  is  less  than  seen  in  the  Space  Shuttle  restorations,  for  two  reasons.  First,  the  HST 
measurement  frames  have  much  higher  quality  than  the  Space  Shuttle  measurement  frames  even  before 
processing.  As  a  result,  there  is  less  room  for  improvement  before  reaching  the  diffraction  limit  of  the 
telescope.  Second,  far  fewer  frames  were  used  for  the  HST  restorations  than  for  the  Space  Shuttle 
restorations. 

6.  CONCLUSIONS  AND  FUTURE  WORK 

We  have  presented  our  blind  deconvolution  algorithm,  called  PCID,  described  its  features,  and 
characterized  its  performance.  Although  many  blind  deconvolution  algorithms  have  been  developed  and 
published  in  the  past,  several  features  of  our  algorithm  extend  the  state-of-the-art  in  blind  deconvolution. 
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First,  we  have  parallelized  the  algorithm  so  that  it  can  run  on  commodity  clusters  -  specialized  hardware 
is  not  required.  We  demonstrated  the  speed  improvements  using  a  Cray  XD-I  cluster  and  showed  how 
PCID’s  scalability  properties  depend  on  the  size  of  the  data  arrays.  We  have  also  developed  Cramer-Rao 
lower  bound  expressions  for  blind  deconvolution  and  carried  out  an  extensive  comparison  study  that 
showed  that  PCID  closely  approaches  these  theoretical  limits.  Finally,  we  presented  image  restorations  of 
the  Space  Shuttle  and  the  Hubble  Space  Telescope  from  data  collected  by  two  separate  telescopes,  with 
and  without  adaptive  optics,  and  showed  that  significant  improvements  in  image  quality  were  obtained. 

Future  work  includes  implementing  the  complete  MACOPT  [43]  conjugate  gradient  minimizer 
which  should  noticeably  increase  PCID’s  execution  speed.  For  further  speed  improvements,  we  plan  to 
implement  the  distributed-memory  version  of  FFTW  when  it  becomes  available  should  it  prove  to  be 
faster  than  our  current  2-D  FFT.  We  will  further  increase  the  automatic  parameter  calculation  properties 
of  PCID  to  make  it  more  user  friendly.  We  are  also  in  the  process  of  writing  a  web-based  GUI  that  will 
be  a  front  end  to  PCID.  In  addition,  this  GUI  will  be  used  for  data  I/O,  job  submission,  and  data 
visualization. 
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